Dynamics of a tracer granular particle as a non-equilibrium Markov process. 
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The dynamics of a tracer particle in a stationary driven granular gas is investigated. We show 
£1h how to transform the linear Boltzmann equation describing the dynamics of the tracer into a master 

equation for a continuous Markov process. The transition rates depend upon the stationary velocity 
distribution of the gas. When the gas has a Gaussian velocity probability distribution function 
. (pdf), the stationary velocity pdf of the tracer is Gaussian with a lower temperature and satisfies 

detailed balance for any value of the restitution coefficient a. As soon as the velocity pdf of the 
gas departs from the Gaussian form, detailed balance is violated. This non-equilibrium state can 
be characterized in terms of a Lebowitz-Spohn action functional W(t) defined over trajectories of 
time duration r. We discuss the properties of this functional and of a similar functional W(t) 
which differs from the first for a term which is non-extensive in time. On the one hand we show 
that in numerical experiments, i.e. at finite times t, the two functionals have different fluctuations 
and W always satisfies an Evans-Searles-like symmetry. On the other hand we cannot observe the 
verification of the Lebowitz-Spohn-Gallavotti-Cohen (LS-GC) relation, which is expected for W(t) 
at very large times r. We give an argument for the possible failure of the LS-GC relation in this 
situation. We also suggest practical recipes for measuring W(t) and W(t) in experiments. 
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I. INTRODUCTION 



. When a collection of macroscopic grains (diameter ranging from 10 cm up to 10 cm) is vigorously shaken, 
' under appropriate conditions of packing fraction, amplitude and frequency of the vibration, a stationary gaseous 
state, usually referred to as "granular gas", can be obtained 0. Kinetic energy is dissipated into heat during 
0^ ■ collisions among grains, and this energetic loss is balanced by the external driving, i.e. the shaking of the container. 
In recent years the interest on granular gases has strongly increased: it has been shown that the whole machinery 
of kinetic theory can be applied on simplified but realistic models, leading to theoretical results in good agreement 
with experiments 0, Q • At the same time a growingly rich phenomenology has emerged from the laboratory, showing 
that a granular gas is a sort of Pandora's box for non equilibrium statistical mechanics. Just to mention the most 
famous and well established aspects of these systems, we recall that a granular gas displays non-Gaussian velocity 
behavior, the breakdown of energy equipartition and of thermodynamic equilibrium, spontaneous breaking of many 
symmetries (shear instability, cluster instability, convection instability, etc.) and a general strong tendency to reduce 
the range of scales available to a hydrodynamic description. Nevertheless, in a carefully controlled environment, one 
can avoid all spatial effects, obtaining a homogeneous granular gas whose main feature is to be stationary and out of 
thermodynamic equilibrium, traversed by an energy current that flows from the external driving into an irreversible 
sink constituted by the inelastic collisions among the grains. In this case it is straightforward to measure a so- 



called granular temperature T g — (|v| /d) (where d is the space dimension) which is different, in general, from the 
temperature T characterizing the external driving. We stress that, in real experiments, it is hard to obtain a stochastic 
way of pumping energy into a granular gas: the simplest situation is a vertically vibrating mechanism, where density 
gradients arise. In this paper our interest goes to models where the action of the external driving is in fact that of 
a homogeneous thermostat. Such a kind of homogeneous heating has already been achieved experimentally with a 
granular mono- layer vibrated on a rough plate Very recently granular gases have been used to probe the latest 
theoretical results concerning fluctuations in non-equilibrium systems. In particular the so-called Gallavotti-Cohen 
Fluctuation Relation (GCFR) [j| has been put under scrutiny 0, Q ■ This relation is a constraint in the probability 
distribution function (pdf) of the fluctuations of entropy flux in the system, which is rigorously demonstrated under 
certain hypothesis when the entropy flux is measured by the phase space contraction rate. Of course an all-purpose 
definition of "entropy flux" in systems out of equilibrium does not exist. The authors of those recent studies have 
tried to use the rate of energy injection from the external driving, i.e. the injected power, as a entropy flux. Even 
if in a verification of the GCFR was claimed, other studies have shown that the situation is more complex and 
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that, in general, the injected power in a granular gas cannot satisfy such relation [8|. The main difference between 
a granular gas and the prototypical system that should obey the GCFR is the invariance under microscopic time 
reversal, which is violated by inelastic collisions. 

In this paper we have tried to follow a simplified line of reasoning, obtaining a recipe to measure a "flux" which 
is constructed ad hoc with the aim of satisfying the GCFR, or at least its Markovian counterpart, which has been 
put forward by Lebowitz and Spohn In the following we will refer to this relation as to the GC-LS Fluctuation 
Relation. Our idea is to bypass the problem of "strong irreversibility" posed by inelastic collisions (i.e. the fact that 
the probability of observing the time reversibility of an inelastic collision is strictly zero), focusing on the evolution of 
a tracer particle Indeed, in a dilute gas, a tracer particle performs a continuous Markov process characterized 

by transition rates which are always defined, i.e. for any observable transition the time-reversed transition has a 
non-zero probability. Therefore, the projection of a TV-body system onto a 1-body system somehow increases its 
degree of time reversibility. 

Our original aim was to find a quantity, for granular gases, which, by construction, verifies the GC-LS Fluctuation 
Relation. Besides the use of a single particle functional is still interesting because it can easily be reproduced in 
real experiments. Moreover, the fluctuations of this functional are much stronger than the fluctuations of some other 
observable which is averaged over a large number of particles, and this is an evident advantage in an experimental 
verification of a Fluctuation Relation. Anyway, we will show by numerical simulations that this verification poses 
unexpected problems, even after having followed carefully the recipe given by Lebowitz and Spohn. Only at very large 
integration times and for situations very far from equilibrium, the GC-LS Fluctuation Relation is near to be satisfied. 
At the same time an alternative functional can be defined, differing from the first one by an apparently small term. 
This second functional always satisfies a relation which we consider the Markovian counterpart of the Evans-Searles 
(ES) Fluctuation Relation [Llll3|. In summary, our investigation has led to a concrete example, ready to become a 
real experiment, where the difference between two Fluctuation Relations, their meaning and their limits can directly 
be probed. 

The structure of the paper is as follows: in section II we define our model and its description in terms of a continuous 
Markov process, i.e. giving the transition rates that enter its Master Equation. In this section we also show that, 
whenever the surrounding gas has a non-Gaussian velocity distribution, the tracer particle violates detailed balance 
and is therefore out of equilibrium. In section III we define the action functionals that are expected to verify the 
GC-LS and ES Fluctuation Relations. In section IV we discuss the results of numerical simulations, showing a way to 
optimally measure the transition rates and finally discussing the measure of the fluctuations of the action functionals 
and their verification of the Fluctuation Relations. We will draw our conclusions in section V. The derivation of the 
Master Equation from the linear Boltzmann Equation, is given in the appendix. 



II. THE MASTER EQUATION. 



We consider the dynamics of a tracer granular particle in a homogeneous and dilute gas of grains which is driven 
by an unspecified energy source. The first essential feature of the gas is its spatial and temporal homogeneity. The 
tracer particle collides sequentially with particles of the gas coming from the same "population" , independently of 
the position and time of the collision. The gas is characterized by its velocity probability density function P(v), 
which in turn is determined by the unspecified details of the model, such as the properties of the driving or of the 
grains. It is well known that the velocity probability density function of a driven granular gas is non-Gaussian: in a 
homogeneous setup the slight departure from Gaussianity is well reproduced by a first Sonine correction (which will 
be made explicit in the following) . This is expected theoretically and well verified in experiments |14| . 

The second essential property required is the diluteness of the gas: it guarantees that the evolution as well as the 
stationary regime of the velocity probability density function P* (v) of the tracer is governed by a Linear Boltzmann 
equation. In this section we give expressions for the tracer transition rates in generic dimension d and for a generic 
inter-particle (short range) potential. This potential is parametrized by a parameter 7 which is the exponent of the 
term |(v — V) • o>| 7 representing the collision kernel in the Linear Boltzmann equation (v and V are the colliding 
velocities while cl> is the unitary vector joining the particles). For example, a value of 7 = 1 corresponds to the Hard 
Sphere potential (which will be studied in the rest of the article, without losing generality in the results), a value of 
7 = corresponds to the Maxwell Molecules case and a value of 7 = 2 to the so-called Very Hard particles. 

The inelastic collisions with the gas particles, which solely determine the instantaneous changes of the velocity of 
the tracer, are described by the simplest and most used inelastic collision rule: 

v' = v-ml±^-[(v-V)-£]u> (1) 
m + M 



where v and v' are the velocities of the tracer before and after the collision respectively, V is the velocity of the gas 
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particle, m and M are the masses of the tracer and of the gas particle respectively, while cl> is the unitary vector 
joining the centers of the two particles. In the following we will assume m = M, but this will not change the generality 
of our results, in view of the following relation 

1 + a 1 + a' , , in — M + 2ma . . 

m — = — - — , with a = — (2) 

m + M 2 m + M w 

The 7-Boltzmann equation for the tracer, in generic dimension d, reads: 



dP,(v,t) v}" 1 



_ 

dt I 



— y" dvi y" dv2 J dw|(vi-v 2 )-w| 7 P»(vi)P(v 2 ) |<5 - V! + i-^[(vi - v 2 ) • w]o>^ - <5(v - vi^ 



(3) 

where the primed integral denotes that the integration is performed on all angles that satisfy (vi — v 2 ) • <2> > 0. The 
mean free path Iq appears in front of the collision integrals. In the 7 ^ 1 cases, to respect dimensionality, a factor 
«0 _7 is a l so present, where vq is the gas thermal velocity. In the following (when not stated differently) we will put 
Iq = 1 and vq — 1, which can be always obtained by a rescaling of time. 

From the analysis of equation © given in Appendix, we find that the evolution of the velocity probability density 
function of the tracer is governed by the following Master equation: 



dP,(v,t) 
dt 



J dv 1 P.(vi)A'(v 1 ,v)- J dviP»(v)A-(v,vi). (4) 



where P*(v) is the velocity pdf of the test particle. The transition rate K(v, v') of jumping from v to v' is given by 
the following formula: 

K(vy) = (t^)^ ^7-l Av l 7 " d+1 / dv 2rP[v 2 (v,v',v 2r )], (5) 

where Av = v' — v denotes the change of velocity of the test particle after a collision. The vectorial function v 2 is 
defined as 

v 2 (v, v', v 2r ) = u 2ct (v, v')<t(v, v') + v 2r , (6) 

where <r(v, v') is the unitary vector parallel to Av, while v 2t is entirely contained in the (d — l)-dimensional space 
perpendicular to Av (i.e. v 2r • Av = 0). This implies that the integral in expression (J5J is (d — l)-dimensional. 
Finally, to fully determine the transition rate (J5J, the expression of t> 2(T is needed: 

u 2ff (v,v') = — 2 — |Av| +v& . (7) 
1 + a 

Note that equation Q is a Master equation for the evolution of a probability density function. We simplify the 
terminology using the name "transition rate" for the function K(v, v'), which actually is a "transition rate density". 
This density has its natural definition in the following limit: 

*(v,vO=lim ^uefi^vO) 

|dv'|->0 |dv'| 

where /C(v -^ue Bdv'(v')) is the probability that the tracer, after a collision, has a velocity u contained in a sphere 
Sdv'(v') of radius dv' centered in the vector v', having a velocity v before the collision. 

A. Examples: Gaussian and first Sonine correction 

If P(v) = ( 2 irT) d / 2 ex P ( — It)' ^ then immediately follows that the transition rate K(v,v r ) reads: 
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In kinetic theory, one of the most used corrections to the Gaussian is the first non-zero Sonine polynomial ap- 
proximation [Til [T5| . This means assuming that P(v) = ^ry/ 2 ex P \~§r)0- + a ^2 ( l|2 /2P)) with S 2 {x) — 



2 x 2 



d+2, 



d(d+2) 



The calculation of the integral needed to have an explicit expression of the transition 



rate is straightforward: 



This leads to 



K(w,w') = 



dv 2T P(v 2 ) 



(1 



7+1 1— 7 



lo 



\Av 



■y-d+1 



a 2 St\vlJ2T)) 



V2ttT 



(10) 



(11) 



B. Detailed balance 



Here, we obtain a simple expression for the ratio between K(v,v') and if(v',v). When exchanging v with v' the 
unitary vector <r change of sign. Furthermore one has that «2o-(v, v') 7^ v 2l7 (v' , v). Finally it must be recognized that 



K(v',v)dv ~~ K(\', v) ' 



(12) 



This equivalence is naturally achieved by considering the definition (|SJ and taking uniformly the two limits |dv| — ► 
and |dv'| — > 0. From all these considerations and from equation JSJ one obtains immediately: 



K{v,v') _ J dv 2T P[v 2 (v,v')] _ P[v 2a (v,v')} 
K(v',v) /dv 2T P[v 2 (v',v)] " P[v 2<7 (v',v)} 



(13) 



We note that this ratio depends only on the choice of the pdf of the gas, P, and not on the other parameters (such as 
7 or a). However in realistic situations (experiments or Molecular Dynamics simulations) P is not a free parameter 
but is determined by the choice of the setup (e.g. external driving, material details, geometry of the container, etc.). 
Introducing the short-hand notation v 2a — u 2(T (v, v'), v' 2a = u 2cr (v', v) and v„ = • <x, we also note that 



+ (v„ + v' a ) 2 - 2v 2rJ (v a + v' a ) , 



(14) 



from which it follows that 



where A 



K) 



' ^2 



(V2af ~ (V' 2a ) 



-A - 2- 



1 



-A 



1 



-A 



(15) 



.,'12 



i.e. the kinetic energy lost by the test-particle during one collision. When a = 1 



then A 2 = — A (energy conservation). From the above considerations it follows that 
• in the Gaussian case, it is found 



^v0 = A 1-aA 
B i4T(v',v) 2T l + a2T 



• in the First Sonine Correction case, it is found 



1 + a 2T 



(16) 



log 



if(v,v') _ 3 



K(y',y) l + a2T 



log- 



jl + a,^ 1 






2T 


} 


jl + a 2 S$ =1 






2T 


} 



(17) 
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In the case where P(v) is a Gaussian with temperature T, it is immediate to observe that 

P*(v)tf(v,v') = P*(v')#(v',v) (18) 

if P* is equal to a Gaussian with temperature T' = j^T < T. This means that there is a Gaussian stationary solution 
of equation 10} (in the Gaussian-bulk case), which satisfies detailed balance. The fact that such a Gaussian with a 
different temperature T' is an exact stationary solution was known from |10| . It thus turns out that detailed balance 
is satisfied, even out of thermal equilibrium. Of course this is an artifact of such a model: it is highly unrealistic that 
a granular gas yields a Gaussian velocity pdf. As soon as the gas velocity pdf P(v) ceases to be Gaussian, detailed 
balance is violated, i.e. the stationary process performed by the tracer particle is no more in equilibrium within the 
thermostatting gas. We will see in section ITTTl how to characterize this departure from equilibrium. 



C. Collision rates 



The velocity dependent collision rate is defined as 



r(v) = J dv'K(v, v') = ^- J dv'da8[(v - v') • a](v - v') • <rP(v') (19) 

where the last passage is true for hard spheres. In the following, for simplicity and for direct comparison with 
numerical results, we will consider only the model of a two-dimensional hard sphere gas. In this particular case, the 
collision rate reads: 

• for a Gaussian bulk: 



r(v) 



W2T 



e 4T 



2 2 * 

(2T + v 2 )I (^-)+v 2 I 1 (^-) 



AT 



'AT' 



(20) 



where I n (x) is the n-th modified Besscl function of the first kind. Note that from equation (|20|l one obtains the 
total collision frequency: 



1 

V 



dvve 2t' r(v) 



lo 



2tt 



-T 



(21) 



which, using Iq = for the mean free path gives uj c — 2V irTna 2 in the case a = 1, i.e. the known result from 
kinetic theory for the elastic hard disks. Finally note that r(v) for the Gaussian case does not depend upon the 
restitution coefficient (while the transition rates and the total collision frequency do). 

• for a bulk with a Gaussian distribution plus the first Sonine approximation: 



r(v) 



l 8V2T 5 / 2 



e"*r (8(1 + a 2 )T 2 - 8a 2 Tv 2 + a 2 v 4 ) 



(2T + v 2 ) I 



AT 



v 2 h 



at) 



(22) 



In this case the expression for lo c is more involved. 



III. NON-EQUILIBRIUM CHARACTERIZATION 



From the previous section we have learnt that the dynamics of the velocity of a tracer particle immersed in a 
granular gas is equivalent to a Markov process with well defined transition rates. This means that the velocity of 
the tracer particle stays in a state v for a random time t > distributed with the law r{v)e~ rtyV ^ t At and then makes 
a transition to a new value v' with a probability r(v)^ 1 K(v, v'). At this point it is interesting to ask about some 
characterization of the non-equilibrium dynamics, i.e. of the violation of detailed balance, which we know to happen 
whenever the surrounding granular gas has a non-Gaussian distribution of velocity. 

To this extent, we define two different action functionals, following @: 
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n(t) 

>°if(v^v 4 ) 



^) = E lQ g S"I^ (23a) 



W W = log ^-7-; — r + > log — r = log — ^— 23b 

where i is the index of collision suffered by the tagged particle, Vj is the velocity of the particle before the i-th collision, 
v[ is its post-collisional velocity, n(t) is the total number of collisions in the trajectory from time up to time t, and 
K is the transition rate of the jump due to the collision. Finally, we have used the notation V(y\ — ► V2 — > ... — > v„) 
to identify the probability of observing the trajectory vi — > V2 — * ... — * v n . The quantities and W^(t) are 

different for each different trajectory (i.e. sequence of jumps) of the tagged particle. The first term log P P ^7 1 ' 1 , in 
' ' ' ^*l v »(t)J 

the definition of W(t), eq. (|23b|l . will be called in the following "boundary term" to stress its non-extensivity in time. 

The two above functionals have the following properties: 

• W(t) = if there is exact symmetry, i.e. if K(y^ — > v i+1 ) = if(v i+1 — > Vj) (e.g. in the microcanonical 
ensemble); W(t) = if there is detailed balance (e.g. any equilibrium ensemble); 

• we expect that, for t large enough, for almost all the trajectories \im s ^ oa W(s) /s — lirn s _,oo W(s)/s — 
(W(t)/t) = (W(t)/t); here (since the system under investigation is ergodic and stationary) the meaning of 
() is intuitively an average over many independent segments of a single very long trajectory; 

• for large enough t: 1) at equilibrium (W(t)) = (W(t)) = 0; 2) out of equilibrium (i.e. if detailed balance is 
not satisfied) those two averages are positive; we use those equivalent averages, at large t, to characterize the 
distance from equilibrium of the stationary system; 

• if S(t) = — J dv P*(v,t) log P*(v,t) is the entropy associated to the pdf of the velocity of the tagged particle 
P*(v, t) at time t (e.g. —H where H is the Boltzmann-H function), then 



±S(t)=R(t)-A(t) (24) 

where R(t) is always non-negative, A(t) is linear with respect to P* and, finally, (W(t)) = J Q dt'A(t'). This leads 
to consider W(t) equivalent to the contribution of a single trajectory to the total entropy flux. In a stationary 
state A(t) = R(t) and therefore the flux is equivalent to the production; this property has been recognized in 

• FRw (Lebowitz-Spohn-Gallavotti-Cohen fluctuation relation): w(w) — 7r(— w) = w where n(w) = 
lim^oo | log f%y (tw) and fw( x ) I s t ne probability density function of finding W(t) = x at time t; at equi- 
librium the FRw has no content; note that in principle n'(w,t) — 7 log f w (tw) 7^ tt(w) at any finite time; 
a generic derivation of this property has been obtained in [9|, while a rigorous proof with more restrictive 
hypothesis is in 16:j; the discussion for the case of a Langevin equation is in [l7| . 

• FRyy (Evans-Searles fluctuation relation): W(w,t) — w(—w,t) = w where W(w,t) — - log f£p(tw) and f^(x) is 
the probability density function of finding W(t) = x at time t; at equilibrium the FR^r has no content; 

On the one side we have called FRw a "Lebowitz-Spohn-Gallavotti-Cohen" fluctuation relation, following |9|, where 
the analogy with the original Gallavotti-Cohen Fluctuation Relation has been stated explicitly. On the other side 
we have called the FR^ a "Evans-Searles" Fluctuation Relation, inspired by the following analogy. The original 
Gallavotti-Cohen Fluctuation Theorem concerns the fluctuation of a functional fl(t) which is the time-averaged phase 
space contraction rate, i.e. 

tfl(t) = [ A(T(s))ds (25) 

where A(r(s)) is the phase space contraction rate at the point T(s) in the phase space visited by the system at time 
s. On the other side, the original Evans-Searles Fluctuation Theorem concerns the fluctuation of a functional 

f2(f) (often called "dissipation function") defined as 
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where /(T) is the phase space density function at time 0, p(SVr(T)) is the probability of observing at time an 
infinitesimal phase space volume of size 5Vp around the point T, and It is the time reversal operator which typically 
leaves the position unaltered and changes the sign of the velocities. The last equality in equation (|2*f)|l follows from the 
Liouville equation written in Lagrangian form: d ^'^ = — A(r)/(I\ t) (details on its derivation can be found in |13|). 
The ratio of probabilities appearing at the end of the r.h.s. of i|26[) can be thought as the deterministic equivalent of 
the ratio of probabilities appearing at the end of the r.h.s of definition l|23b|) . The term log jt^tcj] is instead analogous 
to the term log jr-j^P^ ) , since the system can be prepared at time in such a way that its phase space density 

function / is arbitrarily near to its stationary invariant measure. The analogy between those two terms, and therefore 
between the two pairs of functionals SI, f2 and W, W, cannot be derived mathematically rigorously, at the moment, 
but formally at least, it looks rather convincing. Note also that the difference between the two functionals of each 
couple vanishes when the invariant phase space distribution function becomes constant, e.g. in the microcanonical 
ensemble. 

As a final remark of this section, consider the more general case where the tracer particle feels the external driving 
too: this means that other transitions are possible. The simplest case corresponds to an external driving modeled as 
a thermostat which acts on the grains in the form of independent random kicks |l4j . It is straightforward to recognize 
that the transition rates competing to the velocity changes due to these kicks are symmetric and therefore do not 
contribute to the action functionals in eq. (|23[l . On the other hand, when one wants to include the effect of walls 
(both still or vibrating), new transition rates should be taken into account. A collision with a wall can be modeled as 
a collision with a flat body with infinite mass, so that 

v' a = ~a w v a + (1 + a w )v w (t) (27) 

where a is the direction perpendicular to the wall and v w (t) is the wall velocity in that direction at the time t of the 
collision, while a w is the restitution coefficient which describes the energy dissipation during such collision. If the wall 
is still (v w = 0) and a w < 1 then the transition is non-reversible and this results in a divergent contribution to the 
action functionals. If the wall is still but elastic, then again the transition rates are symmetric and do not contribute to 
the action functionals. Finally, if the wall is vibrating the transition rates will depend upon the probability P w {v w = x) 
of finding the vibrating wall at a certain velocity x. 

IV. NUMERICAL SIMULATIONS 

A Direct Simulation Monte Carlo is devised to simulate in dimension d = 2 the dynamics of a tracer particle 
undergoing inelastic collisions with a gas of particles in a stationary state with a given velocity pdf -P(v). The 
simulated model contains three parameters: the restitution coefficient a, the temperature of the gas T (which we take 
as unity) and the coefficient of the first Sonine correction which parametrizes the velocity pdf of the gas, a-i- The 
tracer particle has a "temperature" < T and a velocity pdf that is observed to be well described again by a first 
Sonine correction to a Gaussian, parametrized by a coefficient a^. We also use Iq — 1. This means that the elastic 
mean free time is rf = (it is larger in inelastic cases). 

The quantities W(t) and W(t) are measured along independent (non-overlapping) segments of time-length r 
extracted from a unique trajectory. 

A. Verification of transition rates formula 

In this subsection we verify the correctness of formula 10 for the Gaussian bulk and formula I|I1|) for the non- 
Gaussian (Sonine-corrected) bulk. Since the algorithm used (the DSMC ^^|) is known to very well reproduce the 
Molecular Chaos assumption, which is the unique hypothesis used to obtain those formula, we consider this verification 
a mere check of the consistency of algorithms and calculations. On the other side, while performing this verification, 
we have devised a way of optimizing the measurement of the transition rates, which can be adopted in experiments. 

The brute force method of measuring K in an experiment or a simulation is to observe subsequent collisions for a 
long time t and cumulate them in a 4-dimensional (2d-dimensional) matrix DY t : this means that DY t (v x ,v yi v x ,v' y ) 
contains the number of observed collisions such that the precollisional and postcollisional velocities are contained in 
a cubic region centered in (v x , v y ,v' x ,v' y ) and of volume dv x x dv y x dv' x x dv' y finite but small. At the same time it is 
necessary to measure P*(v x , v y )dv x dv y for the tagged particle. These quantities are related by the following formula: 

DY t (v x ,v v ,v' x ,v') =tP*(v x ,v y )K(-v,V)dv x dv y dv x dv' (28) 
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which can be immediately inverted to obtain K(v, v'). Anyway this recipe may require a very large statistics, because 
of the high dimensionality of the histogram DY t . 

Then we notice that K is a function of only v a ,v a i. It becomes tempting, therefore, to reduce the number of 
variables from 2d to 2 in order to optimize the procedure. A histogram DY t of dimensionality 2 (in any dimension) 
can be obtained, so that each element DY t (u a , u' a ) contains the number of observed collisions such that the projection 
along the direction a of the precollisional velocity is u a and that of the postcollisional velocity is u' a . We note that 



DY t (u a ,u' a ) = du a du' a J DYt(v x ,v y ,v' x ,Vy)5(u a - v • <t(v, v'))6{u' a - v' • <t(v,v')) = 

du a du' a t J dv x dv y dv' x dv' y K(-v,V)Pf(-v)d(u a — v ■ <r(v, v'))(5(w' CT — v' • <r(v, v')) = 

duadu'^t / dv x dv y / dv' a dv' a J(v x , v y> v a , u^)if(v, v')P*(v)5(u CT - v^S^ - v' a ). (29) 



In the last passage, we have done the change of variables (at fixed v x , v y ) v' x , v' y — > v a (v x , v y ), v' a (v x , v y ). This implies 
the appearing of the associated Jacobian. Moreover, it happens that v-er(v, v') = v a and v' -ct(v, v') = v' a . Postponing 
the problem of finding the Jacobian J , we can absorb the Dirac deltas, obtaining: 

DY t {u a ,v! a )=du a dv! a tK{u a ,v! a ) dv x dv y J(v x ,v y ,u cr ,u' a )P t ,(-v) (30) 



where K(u a , u' a ) is just the transition probability K as a function of u a , u' a (which is its natural representation). The 
change of variables is given by the following rule: 



Ja — v x 



\j K - V X ) 2 + (V' y - Vy) 2 ^ K ~ V X ) 2 + (V' y ~ Vy) 



<^ (31a) 



\j (V x ~ V X ) 2 + (V' y - Vy) 2 V {V' x - V x ) 2 + {V' y - Vy) 

After calculations, the Jacobian reads: 

\ V '<J- V °\ af.fi , „,2 „,2 



(31b) 



j = j g u Q{< +v z_ <y (32) 



v% + vl - v 2 



From the expression of the Jacobian, and taking into account the isotropy of P*(v)dv = uP*(v)dvd0 one gets 

DYiiu^u'g.) =duvdu' ff tK(u a ,u' a .)\u' <r — u ( r\£(ucr), (33) 

where 

1 f f P fvl 

£{u a ) = — - / dv x dv y J(v x ,v y ,u a ,u' cr )P*(v) = / dv x dv y —= * :9{vl + v 2 - u 2 a ) = 



.2 _i_ t,2 — /j 1 2 



2tt / dw = 2tt / d»P.(v^Tuf). (34) 

I Vv 2 -ul 



I] 



As already discussed, the result by Martin and Piasecki [Hj explains that when the velocity-pdf of the gas is 
Gaussian, then also the tagged particle has a Gaussian pdf. In particular if the bulk has a temperature T, the tagged 
particle has a temperature T 1 = §-^jP Then it is easy to find that 



\l + aj ^/2^T 



a J V2^T 

and 
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so that the theoretical expectation for the measured array DY t is 

2 



DY t {u a , u' a ) = du a du' a t 



which for the case a = 1 reads 



\u a - u a \ exp 



1 1 



2T(l + a) 2 



(ul + «) 2 + (a- l)u a u' a ) 



DY t {u a ,u' a ) = du a du' cr t- — \u' a . -u CT |exp 



2T 



(37) 



(38) 



It should be noted that, for any value of a, DY t is symmetric, i.e. DY t {u a , u' a ) = DY t {u' a , u„) and this is the reflection 
of the fact that in this Gaussian case the detailed balance is always satisfied. 

In left and center frames of figure the results of DSMC simulations with Gaussian bulk and a = 1 and a = 0.5 
are shown for some sections of the measured DY t , showing the perfect agreement with the theoretical fits. 




FIG. 1: Measure of transition rates in Monte Carlo simulations. Sections of the function 



for the values u a = —0.1, 
1 and a = 1 



tdu<jdu' a 

Ua = 0.9 and u a = 3.5 for three different cases. Left: the gas has Gaussian velocity pdf with temperature T 
(thermodynamic equilibrium). Center: the gas has a Gaussian velocity pdf and a = 0.5 (detailed balance in absence of 
thermodynamic equilibrium, T„ = 0.6 = T'< T). Right: the gas has non-Gaussian velocity pdf characterized by a first Sonine 
correction with 02 = 0.1 (non-equilibrium stationary state, T, = 0.62 ~ T 1 < T and a| = 0.057 <C 02). The dashed lines mark 
the theoretical fits given in equation 1381 . 1371 and 1391 . 



Long calculations lead to a similar result for the non-Gaussian case, when the gas has a velocity distribution given 
by a Gaussian, with temperature T, corrected by the second Sonine polynomial with a coefficient 0,2- In this case an 
exact theory for the stationary state of the tracer particle does not exist. Nevertheless, numerical simulations show 
that the velocity distribution of the tracer, P*(«) is well approximated by a Gaussian with temperature corrected 
by the second Sonine polynomial with a coefficient a^. In particular we have evidence that T* w V = f±±T (the 
Martin-Piasecki temperature, which is exact in the Gaussian case) and < &2- For hard spheres we obtain: 



DY t (u a ,u' a ) 
du n du't 



1 



1 + aJ 2V7X 



\u a - u a \ exp 



( (u'„ — U a ) 



2T» 



3 3ui 
a 9 1 - 

4T. 



1 + a 2 S: 



d=l 
2 



(tT^Kt -u a ) + u a ^ 
2T 



(39) 



this formula is very well verified by the transition rates observed in simulations, as shown in the rightmost frame 
of figure n 



B. Fluctuation Relations 



In this paragraph we want to show the numerical results obtained measuring W(t) and W(t) in DSMC simula- 
tions |18| . In particular our main aim is to verify the Fluctuation Relations FRw and FR^. Actually there are three 
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equivalent relations that can be written at finite times r. 



as T (W) 



log r{w)- log r(w)=w 

-[log/ T (rw) - log/ T (-ru))] = w 

T 

— l — [log f T (T(w)q) - log f T (-T(w)q)} = q 



(40a) 



as*(w) 



(40b) 



as**(q) 



(40c) 



where f T (x) is the probability density function of finding one of the two functional after a time r equal to x (i.e. 
or defined at the beginning of this section), while (w) is a long time average of W(t)/t. The verification of the 
GC-LS Fluctuation Relation requires the measure of as**(q) at very large values of t, for values of q at least of order 
1, which corresponds to measure as T (W) for values of W at least of order t(w) at large values of r. On the contrary 
the verification of the ES Fluctuation Relation requires equation l|4()a|l to be true at all times. We have chosen to 
display as T (W) vs. W as well as as T (W) vs. W. 

To measure the quantity W(t) on each segment of trajectory it is assumed that P*(V) is a Gaussian with a 
first Sonine correction and the values T* and a?, are measured during the simulation itself and used to compute the 
"boundary term" appearing in the definition of W, eq. I|23b|) . All the results are shown in figure [21 for many different 
choices of the time r and of the parameters a and 02- We display in each figure the value of (w) = lim^oo W(t)/t. 
Following the idea that this number is an average entropy production rate, we use it to measure how far from 
equilibrium our system is. This idea is well supported by the results of our simulations: (w) is zero when 02 = and 
increases as a is decreased and is increased. We recall that lmv^oo W(t)/t = lim T _ ) . 00 W(t)/t = (w), since the 
difference between the two functionals has zero average at large times r. We also remark that the distribution of both 
quantities W and W are symmetric at equilibrium (i.e. when 02 = 0). 

We first discuss the results concerning the fluctuations of W(t), identified in figure [3 by squared symbols. The pdf 
of these fluctuations are shown in the insets. They are strongly non-Gaussian, with almost exponential tails, at low 
values of r for any choice of the parameters. At large values of r (many hundreds mean free times), the situation 
changes with how far from equilibrium the system is. At low values of (w) the tails of the distribution are very similar 
to the ones observed at small times. At higher values of (w) the distribution changes with time and tends to become 
more and more Gaussian. We must put a great care in the interpretation of this observation. There is not a firm 
guiding principle to determine how large the time t and the values of W should be, in order to claim that "large 
deviations" are being probed. A simplistic approach is the following: every time one observes non-Gaussian tails then 
one concludes that these tails are, in fact, large deviations, since "small deviations" should be Gaussian. Anyway this 
approach is misleading: a sum of random variables can give non-Gaussian tails if the number of summed variables is 
low and the variables itself have a non-Gaussian distribution. Therefore time plays a fundamental role and cannot be 
disregarded. The observation of non-Gaussian tails alone is useless. We will discuss in the next section the problems 
related to the true asymptotics of the pdf of W(t). Here we stick to the mere numerical observations: we have not 
been able to measure negative deviations larger than the ones shown in figure |21 Distributions obtained at higher 
values of r yielded a smaller statistics and very few negative events. Using the distributions at hand, we could not 
verify the GC-LS Fluctuation Relation. In the last case, which is very far from equilibrium (02 = 0.3, a = 0), a 
behavior compatible with the GC-LS Fluctuation Relation is observed, i.e. as T {W) « W, but the range of available 
values of W is much smaller than (W). At this stage, and in practice, we consider such results a failure of the GC-LS 
Fluctuation Relation for continuous Markov processes. 

On the other hand, things for the functional W(r) are much easier. The pdfs of W(t) always yield deviations from 
the Gaussian, but weaker than those observed for W(t). The ES Fluctuation Relation is always fairly satisfied in all 
non-equilibrium cases at all times r. The verification of both Fluctuation Relations is meaningless in the equilibrium 
cases (0,2 = 0) because as T = 0. 



As already noted, the two action functionals defined in formula ill'MI) differ for the presence/ absence of a boundary 
term. In principle one should expect that it is not possible to distinguish between the large deviations (i.e. the leading 
factor of the pdf at large times) of W(t) and W(t). Anyway there are no arguments to predict the threshold time 
above which the two functionals (i.e. most of their fluctuations) become indistinguishable from the point of view of 
large deviations Moreover, to our knowledge, it is not proved that the large deviation functions of W and 

W coincide. In this paragraph we first discuss numerical observation: they concern, generically, "deviations" , i.e. 
fluctuations at finite time and finite values of the variables. Finally we also argue about the asymptotic behavior of 



C. Difference between the fluctuations of W(t) and of W(t) 
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FIG. 2: Verification of the Fiuctuation Reiations for W(t) (squares) and W(t) (circles). Each line is composed of two graphs 
and shows the results for a particular choice of a and ai (the coefficient of the first Sonine correction characterizing the non- 
Gaussianity of the velocity pdf of the gas): the left graph is at small times and the right graph at large times. In each frame 
the inset contains the pdfs of W and W, while the main plot shows as T (W) vs. W as well as as T {W) vs. W . The times r 
are rescaled with the tracer mean free time r c (which varies with the parameters). For each choice of parameters, inside the 
right graph, it is shown the value of (W(t)) and of (w) = (W(t))/t = (W(t))/t = lim T ^oo W(t)/t = lim T ^oo W(t)/t which 
marks the distance from equilibrium. The dashed line has slope 1. 



12 



the pdfs of W and W . As already stated, in numerical experiments, averages have been always obtained using many 
independent segments of length r from a very long trajectory. This amounts to sampling the stationary velocity pdf 
of the tracer. 




FIG. 3: Cumulants of the fluctuations of W(t) (empty symbols) and W(t) (full symbols). A): average values of W(r) and W(t) 
as a function of r and for different choices of the parameters a and 0,2- B) second cumulants of W{t) and W(t) as a function 
of r. C): ratio between the second cumulant of W(t) and that of W(t) as a function of r. D): ratio between 5y (W(r) 2 ) c 
and (W(t)). When this ratio is much smaller than 1 the probability of observing negative events in the fluctuations of W(t) 
becomes extremely low. All the times r are rescaled by the tracer mean free time. 

In figure [3] we have summed up the results of numerical simulations with different values of a, 02 and r. The 
averages (W(t)) and (W(r)) converge to the same value in a time smaller than 100 average collision times (frame 
'A' of figurc[3J). The analysis of fluctuations (frames 'B' and 'C'), by means of the measure of their variance, instead 
indicates that the convergence is much more slower. Remarkably the convergence is slower when the system is closer 
to equilibrium. We also tried to answer to the question: is it possible to verify the Fluctuation Relation FRw at 
finite times in the regime (large r) where the two functionals have negligible differences? The main obstacle to the 
verification of a Fluctuation Relation is the lack of negative events at large values of r. The probability of finding 
a negative event is equivalent to the probability of finding a negative deviation from the mean of the order of the 
mean itself. We have seen that when the pdf of W(t) seems to converge to its asymptotic scaling form, its tails are 
not far from those of a Gaussian and the bulk (i.e. its "small deviation" range) is clearly Gaussian. This means that 
there is a very sensitive decay of probability in the range of a finite number of standard deviations from the mean. 
From a numerical point of view it is very unlikely to observe values of W(t) below the value (W(t)) — b^J (W 2 (t)) c , 
where (W m (T)) c is the m-th cumulant of W(r), or, equivalcntly, that it is very unlikely to observe negative events if 
5y/ {W 2 {t)) c <C (W(t)). Frame 'D' of figure |3| shows the ratio between 5y/ '(W 2 (t)) c and (W(t)). The conclusion of 
this analysis is that when r is large enough to guarantee, at least at the level of the second cumulant, the convergence 
of the pdf of W (or its convergence to the pdf of W), the probability of observing negative events is so small that a 
huge statistics is required, in order to make reliable measurement of the Fluctuation Relation. 

We now discuss one possible reason for such a slow convergence of the pdf of the fluctuations of W(t). We will call 
it the "boundary term catastrophe". As a matter of fact it seems that W(t), even at very large values of t, more 
evidently in the close-to-equilibrium cases, remembers its own small time fluctuations. This is clearly seen in figureEl 
four topmost frames. The tails of W(t) at large r are almost identical to those at small r. What is happening? If 
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one looks closely at expression (11711 . which is summed over many collisions to give Wit), discovers that 

|2 -' 12 v- ■ 

+ > log — 

Q2- 



l + a 21 ' l + aa — 



where we have used a shortened notation to identify the "non-equilibrium" part which is the sum over all the collisions 
in the time r of the logarithms of ratios between the Sonine contributions. The sum of the energy differences reduces 
to a difference between the first and the last term, which we can call again a "boundary term" . It easy to realize that 
the pdf of this term alone, being the pdf of a difference of energies, has exponential tails (in granular gases they can be 
even slower). This is the dominant term in the fluctuations of W(t) at small r and it can dominate the fluctuations 
of W(t) even at very large r, depending on the amplitude of the fluctuations of the "Sonine" term. The "Sonine" 
term is small near equilibrium and increases as the system gets far from equilibrium. This explains why the memory 
is stronger (i.e. convergence is lower) as the system is closer to equilibrium. 

The situation dramatically changes when looking at the second functional, W(t). The "boundary term" discussed 
above almost perfectly annihilates with the boundary term appearing in the definition of W(t) . Following the numerical 
observation that the pdf P» of the velocity of the tracer is almost Gaussian with temperature T* « T' = jr^T, in 
fact we see that 

P.(vi) ^ -|vi| 2 + K(d 2 „ 3-a-|vi| 2 + |v»( t) | 2 
-P*(v„ (t) ) ~ 2T„ ~l + a 2T 1 ' 

The error in this approximation is of the same order of one of the Sonine terms appearing in the sum (|41|) and is 
therefore subleading at large times. This crucial observation explains the absence of exponential tails in the pdf of 

W(t). 

It should be stressed that this "boundary term catastrophe" has relevant consequences not only in the realm of 
numerical simulations, but also from the point of view of large deviation theory. We suspect that, even with a 
computer of infinite power, i.e. with the ability of measuring with infinite quality the pdf of W(t) at any time r, the 
results of a large deviation analysis would yield similar results: the exponential tails related to the pdf of the v\ — v„ 
are never forgotten. This is the consequence of the following simple observation: following the definition of "large 
deviations rate" of a pdf f t which depends upon a parameter t, we can also calculate the large deviations rate of a 
pdf F which does not depend upon that parameter. We recall the definition: 

n{w) = lim -log/tO). (43) 

t^oc t 

It can be immediately seen that such definition gives a finite result if applied to a pdf F with exponential tails, even 
if it does not depend upon t. This is equivalent to state that exponential tails always contribute to large deviations. 
The apparent loss of memory that can be observed in the strong non-equilibrium cases of figure [21 is only due to 
the limited range available in the measure of the pdf: these tails, in fact, get farther and farther as time increases. 
Moreover we recall that such exponential tails are associated to the distribution of energy fluctuations and are always 
present in a system with a finite number of particles, not only in the 1-particle case. They disappear only in the 
thermodynamic limit or in the presence of physical boundaries. We conclude this discussion mentioning that the 
boundary term catastrophe is responsible for the failure of GC Fluctuation Relation in other systems [Hj, |2(J, l2ll |22| . 



V. CONCLUSIONS 



In conclusion we have studied the dynamics of a tracer particle in a driven granular gas. The tracer performs 
a sequence of collisions changing its velocity, thus performing a Markov process with transition rates that can be 
obtained from the linear Boltzmann equation. We have given a general expression for these rates, and we have also 
calculated them for a number of possible interesting cases. There is a major discrimination between the Gaussian 
and the non-Gaussian case, depending upon the velocity distribution of the gas particles. In the Gaussian case the 
tracer is at equilibrium, even if the collisions are dissipative and energy equipartition is not satisfied. In the non- 
Gaussian case the tracer reaches a statistically stationary non-equilibrium state: in this state a dissipative flow can 
be measured in the form of an action functional that has been defined by Lebowitz and Spohn . The average of 
this flow well describes the distance from equilibrium. Moreover the rate of variation of the entropy of the tracer 
can be split up in two contributions, one always positive and the second given by the average of this dissipative 
flux. Following the Lebowitz and Spohn "recipe" we have tried to verify the Gallavotti-Cohen fluctuation relation for 
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the action functional, realizing that this verification is not accessible. At the same time we have perfectly verified a 
second fluctuation relation, for a slightly different functional, which we have called Evans-Scarlcs Fluctuation Relation 
because of a formal analogy with it. We have discussed the possibility that a non-extensive term (in time) in the 
definition of the Lebowitz-Spohn functional can spoil the applicability of the large deviation principle. In our numerical 
simulations we have also very well verified the correctness of the formula for the transition rates, showing a way to 
optimize their measure in real experiments. A key point of our work is in fact that such Lagrangian approach in 
granular gas theory should be applied to experiments on dilute granular matter, in order to test directly the limits 
of the different Fluctuation Relations. Future work will be directed to the study of transients and non-homogeneous 
situations. The Evans-Searles Fluctuation Relation is expected to hold also in non stationary regimes, such as the 
transient relaxation of a packet of non-interacting tracer granular particles prepared in some non-typical initial state. 
Non-homogeneous situations, such as a boundary driven granular gas, are characterized by the presence of spatially 
directed fluxes. It will be interesting to study the "action functional" defined in this work and compare it with other 
measures of transport. 

Acknowledgments.- The authors acknowledge A. Barrat for several fruitful discussions. A. P. warmly thanks L. 
Rondoni for discussions and a careful reading of the paper, and acknowledges the Marie Curie grant No. MEIF-CT- 
2003-500944. 



Appendix. Derivation from the linear Boltzmann Equation. 



We first discuss in detail what happens in a collision and then give a rigorous derivation of the master equation. 
The collision rule for inelastic hard spheres reads: 



vi = vi — (V12 • ar)cr (44) 

where <r is the direction joining the centers of the two colliding particles. There are some consequences of the collision 
rules which have to be remarked. For simplicity we assume to be in dimension d = 2. 

• Avi = v[ — vi is parallel to <r, i.e. 9 = arctan where a x = cosd and a y = sind. This is equivalent to 
recognize that the velocity vi changes only in the direction a. The fact that vi 2cr must be negative determines 
completely the angle 9, i.e. the unitary vector a. From here on, we call Avi = Avi<j = Avi • a. 

• v 2a = v 2 • d- = T ^Av 1 + v la = j^v' lrT - ±=%v lr7 . 

• From the previous two remarks, it is clear that Av± determines univocally a and v 2l7 . The component of v 2 
which is not determined by Avi is the one orthogonal to a. We call f the direction perpendicular to a, i.e. the 
vector of component (— sin 9, cos 9). We define v 2t = v 2 • t. 

From the above discussion, it is easy to understand that the transition probability for the particle 1 to change 
velocity during a collision, going from Vi to v' x must be 



K(w 1 - vi) = <7(vi,vi) j dv 2T P(v 2 ) (45a) 

v 2 = v 2rT a + v 2t t (45b) 

it = (cos 9, sin 9) (45c) 

f = (-sin (9, cos (9) (45d) 

6 = arctan (45e) 
Avi x 

a || vi - vi (45f) 
2 

Via = — — Avt + v la (45g) 



where P(v) is the 1-particle probability density function for the velocity in the bulk gas. The constant of proportion- 
ality C must be of dimensions 1/ length so that K has dimensions 1 / (velocity d time) which is expected because K is 
a rate of change of the velocity pdf (in d dimensions). 
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Now, we want to obtain the complete result, that is rigorously transform the usual linear 7-Boltzmann equation 
for inelastic models in a Master Equation for a single-particle Markov chain, i.e. as 



dt 



J dviP*( Vl )i<:(vi -> v) - J dviP*(v)if(v -> vi). (46) 



with P*(i>) the velocity pdf of the test particle. 

From this definition (and from the usual linear Boltzmann equation for models containing a term |vi2 • w| 7 , 
equation it follows that 



K (vi v[) = ^ — / dv2 y d ^ e ( v i2 • ^)l v i2 • o>| 7 P(v 2 )<5 jvi - Vi + ^-^[v 12 • a>]wj . 



(47) 



where P(u) is the velocity pdf of the bulk gas. Using that for a generic d-dimensional vector r — rr one has 
6(r — ro) = -j=xS(r — ro)5(r — fo), the previous expression can be rewritten as: 



K{V! - vi) = |dv 2 | du,6(v 12 • h>) l^lgH p( Va )g(g + w)<S f At, + ^|v 12 • , (48) 

where Aw and <r are defined by Vj — Vi = At; cr. Then, performing the angular integration over c2>, one obtains: 

K(yi -> vl) = \^ / dv 2 6(v 12 • a)^-^P(v 2 )«S ^At- + ^|v X2 • erf) . (49) 

Denoting by u 2(T the component of v 2 parallel to er, and by v 2r the (d — l)-dimensional vector in the hyper-plane 
perpendicular to cr, the above equation is rewritten as; 

K(y! - vi) = ^ y dt; 2CT dv 2T e(v 12 • a)^-^P{v 2<J , v 2r )S Uv + ^|v 12 • d>|) . (50) 

Finally, integrating over dv 2(T , Eq. (jSJ) is easily recovered. 
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